Structure–Property Relationship of Piezoelectric Properties in Zeolitic Imidazolate Frameworks: A Computational Study

Metal–organic frameworks (MOFs) are a class of nanoporous crystalline materials with very high structural tunability. They possess a very low dielectric permittivity εr due to their porosity and hence are favorable for piezoelectric energy harvesting. Even though they have huge potential as piezoelectric materials, a detailed analysis and structure–property relationship of the piezoelectric properties in MOFs are lacking so far. This work focuses on a class of cubic non-centrosymmetric MOFs, namely, zeolitic imidazolate frameworks (ZIFs) to rationalize how the variation of different building blocks of the structure, that is, metal node and linker substituents affect the piezoelectric constants. The piezoelectric tensor for the ZIFs is computed from ab initio theoretical methods. From the calculations, we analyze the different contributions to the final piezoelectric constant d14, namely, the clamped ion (e140) and the internal strain (e14int) contributions and the mechanical properties. For the studied ZIFs, even though e14 (e140 + e14int) is similar for all ZIFs, the resultant piezoelectric coefficient d14 calculated from piezoelectric constant e14 and elastic compliance constant s44 varies significantly among the different structures. It is the largest for CdIF-1 (Cd2+ and −CH3 linker substituent). This is mainly due to the higher elasticity or flexibility of the framework. Interestingly, the magnitude of d14 for CdIF-1 is higher than II–VI inorganic piezoelectrics and of a similar magnitude as the quintessential piezoelectric polymer polyvinylidene fluoride.


INTRODUCTION
Energy harvesting is a process where ambient energy in the environment such as mechanical vibrations, light, and heat can be stored and converted into electrical energy. This has been of significant interest to researchers due to its ability to achieve self-powered low-power electronic devices for the Internet of Things and can pave a path to alternate sustainable sources of power. Among the different sources of available energy, kinetic energy in the form of mechanical vibrations is widely available in the environment and various human activities. It can be harvested by piezoelectric materials. The phenomenon of piezoelectricity describes the coupling between the mechanical and electrical properties of materials. This is observed in noncentrosymmetric crystalline structures, where either an electric dipole moment is generated on the application of mechanical stress or a mechanical strain is induced on the application of the electric field. 1 Piezoelectricity has been found in traditional materials such as inorganic ceramic oxides and organics polymers, among which a high piezoelectric constant was found in lead zirconate titanate (PZT, 360 pC/N) and barium titanate (BaTiO 3 , 191 pC/N) for inorganics, while among soft materials, polyvinylidene fluoride (PVDF) and its copolymers (−40 pC/N) have a high piezoelectric constant. Ceramics and polymers have their own set of advantages and disadvantages as piezoelectric materials, which are summarized and reviewed in the previous literature. 2,3 Ceramics have strong piezoelectric properties but are stiff and brittle, whereas polymers have mild processing conditions and excellent flexibility, making them suitable for physically flexible electronics, but do not have very high piezoelectric coefficients. 3,4 Moreover, in contrast to inorganics, the acoustic static impedance of polymers is better matched to that of water and living tissues, making them better equipped to harvest energy in these media. 5 With the aim of combining the best of both polymer and ceramic piezoelectric components, composite and hybrid materials were fabricated. For example, in recent years, hybrid organic−inorganic perovskite materials emerged as potential candidates with high piezoelectric response close to conventional ceramic oxides. One of the initial works on trimethylchloromethyl ammonium trichloromanganese (TMCM-MnCl 3 ) hybrid perovskites (a BaNiO 3 -like structure) shows a piezoresponse d 33 as large as 185 pC/N. 6 These hybrid perovskites have features from both organic and inorganic materials and have the advantages of being easy to process, lightweight, and structurally tunable. Other recent works on hybrid perovskites show the huge potential of these materials in terms of their piezoelectric constants. 7,8 Additionally, due to the regulations imposed on the use of hazardous substances such as lead and other heavy metals in electrical and electronic equipment, there has been a need for exploring lead-free environmentally friendly piezoelectrics. Even among the traditional ceramics, lead-free ceramics such as modified potassium sodium niobate (KNN), bismuth sodium titanate (BNT), and other bulk materials were explored as alternatives to PZT. 9,10 Yet, the disadvantages of ceramics are still relevant for these lead-free materials. Another step toward lead-free piezoelectrics was also achieved by obtaining a large piezoelectric response in metalfree organic perovskites: the d 15 of metal-free MDABCO-NH 4 -X 3 (with X = Cl, Br, or I) calculated from first-principles calculations show large values of 119, 248, and 178 pC/N respectively. 11 Like hybrid perovskites, metal−organic frameworks (MOFs) have a major advantage of structural tunability; the building blocks of MOFs, that is, the metal node and connecting organic linker, can be chosen and varied to obtain lead-free piezoelectric structures.
MOFs are materials that consist of metal ions or inorganic clusters connected by organic linkers through directional coordination bonds into a three-dimensional crystalline nanoporous framework. They are distinguished for their permanent porosity and high surface areas. By proper selection of the metal nodes and organic linkers in MOFs, structural control can be achieved to obtain desired properties for the target application. The pore size of MOFs can be changed from several angstroms to nanometers by varying the length of organic linkers in the MOF. Because of these features, MOFs are ideal candidates for applications including but not limited to gas separation and storage, catalysis, and biomedical imaging. 12,13 In addition to these applications, the secondorder nonlinear optical properties via the rational design of non-centrosymmetric MOFs have also been explored previously. 14−16 So far, only two papers have shown piezoelectric response experimentally in MOFs by quasi-static measurements (measured using a PM200 piezometer): in one d 33 value of 60.10 pC/N was measured for MOFs based on Cd-[imazethapyr]; 17 in the other, the d 22 value of Mn-/Co-based homochiral coordination frameworks is measured to be 6.9 pC/N. 18 There are other works on piezoelectricity in MOFs where electromechanical (EM) response was measured using piezoelectric force microscopy (PFM). 19−21 Whereas they are worthy to be mentioned, it must be noted that several nonpiezoelectric effects can induce additional contributions to PFM response, thus leading to misinterpretation. Several issues due to the calibration of the microscope and challenges in excluding the contribution of electrostatic and other forces have been discussed in the literature, thus indicating that piezoelectric constants with high precision cannot be measured from PFM in most cases. 22−24 Recently, also the theoretical prediction of the piezoelectric tensor for three MOFs whose d ik values range between 2 and 23 pC/N has been reported. 25 Another factor that is considered important in piezoelectric energy harvesting applications is the dielectric constant of the material, as noted from the figure of merit (FOM) with units m 2 /N, where ε 33 T is the relative permittivity (or dielectric constant) of the material at constant stress (T). 26 In relation to the dielectric constant of MOFs, because of their inherent permanent porosity, they generally have low dielectric constants (ε r ) ranging between 1.3 to 6. 27−32 Note that the computationally obtained dielectric constants are static dielectric constants and that the experimental determination via dielectric spectroscopy is often done on air−MOF composites (pellets), both leading to some underestimation of the dielectric constant. Yet, the reported values are very low, and, for example, a study where the determination of the dielectric constant of ZIF-8 was performed via ellipsometry on a continuous film, a method not prone to underestimation, reports a low value of 2.3. 27 These values are very low compared to the dielectric constants of BaTiO 3 , PZT, and PVDF and copolymers shown previously. If a high EM response can be achieved by tuning the topology and building blocks of the material, combined with such low dielectric constants, MOFs present an interesting class of materials to explore for efficient energy harvesters. The experimental and computed data for MOFs are thus quite promising in terms of the piezoresponse. However, the MOFs studied differ widely in their building units and structure. Thus, it is hard to derive a structure−property relationship that will guide the design of better-performing piezoelectric MOFs. In the present work, we aim at filling this gap. To this purpose, we chose to focus on zeolitic imidazolate frameworks (ZIFs), a subfamily of MOFs consisting of M− Im−M, where M is the metal cation (i.e., Zn, Cd, or Co) and Im is the imidazolate organic ligand and its derivatives. The reference to zeolites in the naming of ZIFs relates to their topological similarity. Zeolites consist of tetrahedrally bonded Si 4+ (or Al 3+ ) atoms connected by oxygen atoms, forming Si/ Al−O−Si/Al angles of 145°. Similarly, in ZIFs, the metal (M 2+ ) cations are tetrahedrally coordinated to nitrogen at 1,3positions of the organics imidazolate linkers, leading to a M− Im−M angle of 145°. This leads to this class of metal−organic frameworks, forming the same net topologies as zeolites. Due to the organic linkers and metal−organic bonds, ZIFs generally have higher flexibility compared to zeolites. ZIFs are also known to have excellent thermal (can withstand up to 550°C) and chemical stability. 33 In this work, we specifically focus on sodalite (sod) ZIFs with Zn or Cd metal nodes and Rsubstituted imidazolate linkers (R = CH 3 , Cl, CHO, NO 2 ) as shown in Figure 1. Examined ZIFs belong to noncentrosymmetric space groups (i.e., I4̅ 3m) with a large degree of tunability while keeping the same topology of the prototypical ZIF, that is, ZIF-8 (with Zn and R = CH 3 ). We calculated the pore size, which is the diameter of the largest sphere that will fit into the ZIF framework. The largest pore size for Zn−ZIFs in this work (i.e., ZIF-8, ZIF-90, ZIF-Cl, and ZIF-65) ranges between 11.5 and 12 Å. 33 For Cd−ZIFs (i.e., CdIF-1 and CdIF-8), the pore size varies between 14.2 and 15.2 Å, 34 which is slightly higher than that of Zn−ZIFs. This is due to the slightly longer Cd−N bond compared to Zn−N. Such high porosity should contribute to a lower dielectric constant favorable for piezoelectric energy harvesting. Indeed the total porosity for ZIFs in this work varies between 48 and 60% of the volume of the unit cell. Notably, from an application point-of-view, among the studied ZIFs, ZIF-8 was already shown to have a very low dielectric constant of 2.3− 2.45 with a low loss factor, 27,35 and can be grown into thin films showing good processibility. Overall, the selected ZIFs are ideal candidates for studying structure−property relationships.
In this paper, we present a systematic computational study on how piezoelectric response varies with changes in building units of ZIFs. For this, we change the metal node (Zn, Cd) and substituent on the imidazolate linker in ZIFs. Piezoelectric coefficients for the investigated ZIFs were calculated using different density functional theory (DFT) methods. To understand the contributing factors of piezoelectric response in ZIFs, we dissect the different contributions of piezoresponse to a mechanical strain into clamped ions, internal strain contributions, and Born effective charges (BECs). We then discuss the results of mechanical properties (specifically the compliance constant s 44 ) and the piezoelectric coefficient for the different ZIFs. For this class of materials, we show that the mechanical properties have a dominant contribution to the piezoelectric constant d than the contribution from the piezoresponse e. In the last part of the paper, we compare these results of ZIFs with some existing Zn/Cd-based inorganic and organic piezoelectrics.

THEORY
Piezoelectricity is mathematically described in the IEEE standard for piezoelectricity 36 by a set of four constitutive equations that describe the response of a piezoelectric material to a mechanical load (stress/strain) and electric fields. The relevant set of equations involving the mechanical and electrical variables and their units are discussed in detail in Supporting Information Section 1. The piezoelectric tensors with units [C/m 2 ] and [pC/N], respectively, and of importance in piezoelectric energy harvesting applications to store as energy is the value of d which effectively links the applied stress (T) to the electric displacement (D). We use the Voigt notation for the third-order piezoelectric tensors (e ikl and d kij ) and fourth-order compliance tensor (s ijkl E ) where the indices are given in a compressed matrix notation instead of the tensor notation. Hence, they can be rewritten as e iq , d ip , and s qp E . The piezoelectric constant d ip can be computed from the piezoelectric constant e and the elastic compliance constant s using d ip = e iq s qp E . Computationally, the piezoelectric constant e iq is calculated as the first derivative of the magnitude of the The piezoelectric tensor e can be further separated into two parts: clamped ion (e 14 0 ) and internal strain (e ik int ) contributions. This is according to a well-known scheme used for many inorganic piezoelectrics to understand the origin of piezoelectricity in the material. 37−39 As indicated in the work by Catti et al. 40 for a general case, piezoelectric constant e is shown in eq 1.
Essentially, the clamped ion contribution (e ik 0 ) is the change in polarization P i due to the reorganization of electron density with external strain ε k , while the fractional coordinates x are kept constant in the strained unit cell, that is, ( ) P x i k . The internal strain contribution (e ik int ) provides the change in polarization P i where the atoms are allowed to relax in the strained unit cell in response to the strain. This means that the internal strain contribution and the clamped ion contribution are typically opposite in sign. For the internal strain contribution, e ik int , we sum over all the atoms "s" in the unit cell along three directions "j" (1, 2, and 3). In this contribution, the change in polarization P i with fractional coordinate x sj = 0,

is multiplied by the change in fractional coordinates
x sj with the lattice strain, also termed the relaxation coefficient Essentially, the internal strain contribution refers to the inner deformation of the crystal structure to keep the energy minimum. Moreover, 40 The BEC is a dynamical charge introduced by Max Born and Maria Goppert Mayer. It refers to a change in polarization induced by an atomic displacement under the condition of zero macroscopic electric field. 41 Hence, for a cubic system, such as the investigated ZIFs, the internal strain contribution can be written in terms of BECs, as indicated in eq 2.
where a is the lattice parameter and V is the unit-cell volume of the unstrained structure. Thus, the internal strain contribution is obtained by summing the product of BECs, a second-order tensor (BECs Z s,ij * ) and respective relaxation coefficient ( ) From eqs 1 and 2, fully understanding the e ik in any piezoelectric material requires the determination of three quantities (a) e ik 0 clamped ion contribution, (b) Z s,ij * BECs, and (c) x sj k , the relaxation coefficient, for all atoms in the unit cell.
Note the relevance of the BECs: the high piezoelectric response of the best-performing ceramics is due to very high anomalous dynamical BECs of their transition metals (M) (vide infra). Further research on these dynamical charges in ceramic oxides indicated the mixed covalent and ionic character of the M−O bonds. Dynamical changes in the hybridization of these bonds with atomic displacement results in charge reorganization, which is observed as the BEC of atoms. 41 The piezoelectric tensors for the investigated ZIF space group I4̅ 3m have only one independent piezoelectric constant: e 14 and d 14 . Specifically, for these space groups, eqs 1 and 2 can be rewritten as shown in eq 3. The relation between piezoelectric constants e and d via compliance constant s is shown in eq 4.
where c 44 is the shear elastic constant and s 44 is the compliance constant.

METHODS
All the calculations reported in this paper are performed using the ab initio periodic code CRYSTAL17 42 based on the atom-centered Gaussian basis set. Triple zeta basis sets with polarization functions for all atoms 43 55 were used for structural relaxation. The deviation between experimental and simulated lattice parameters for all functionals is always below ±2%, as shown in Table S1. Using the optimized structure as starting geometry, full piezoelectric and elastic tensors were calculated using the numerical approach based on the geometry optimization of atomic positions at strained lattice configurations. Piezoelectric constant e is evaluated using the Berry phase approach in the modern theory of polarization 56−58 as a first derivative of the Berry phase 59 with respect to strain. 60 This is done in CRYSTAL17 by applying finite strains to the lattice as finite differences of polarization at strained configurations. The elastic or compliance constants c or s are expressed as second derivatives of energy with respect to pairs of strains, wherein CRYSTAL17, the first derivative, is computed analytically, while the second derivative is evaluated as numerical finite differences between strained structures described in the work by Erba et al. 61 The other piezoelectric constant d is evaluated from the relation between e and c or s shown in eq 4. ZIF structures were subjected to a finite strain of ±1.5% for the numerical differentiation of the analytic cell gradients. We obtained piezoelectric and mechanical properties for all five DFT functionals; however, our most elaborate discussion in the main paper will involve the results using B3LYP, which will be motivated later in the text. The raw data files with the representative input files and output files of the calculations done in this work can also be accessed by others at 4TU.ResearchData. 76

RESULTS AND DISCUSSION
In this study, we consider six different ZIFs with Zn/Cd as the metal node and substituted imidazolate as the linker, namely, ZIF-8, CdIF-1, ZIF-90, ZIF-Cl, ZIF-65, and CdIF-8. The four different substituents of imidazolate considered are CH 3 , CHO, Cl, and NO 2 . Figure 2 shows the equilibrium structure of methyl and nitro-substituted ZIFs, where part of the linkers are oriented along the dipolar direction x (denoted 1 in Voigt's notation), the direction along which we model the change in the electric response while applying a shear deformation along the depicted y and z directions (i.e., 2 and 3 directions denoted as 4 in Voigt's notation). This is also the same for other studied ZIFs. We first present the detailed outcomes of e 14 for all ZIFs in terms of the clamped ion (e 14 0 ) and internal strain contributions (e 14 int ). Later, we will discuss the mechanical properties, that is, s 44 and resultant d 14 . The analysis of piezoelectric constant "e" by separating into the clamped ion and internal strain contributions is also applicable to perovskite materials and has been done previously for inorganic piezoelectrics. 37,39,40,62 We used a series of different functionals (B3LYP, PBE0, M06, M06L, and M062X) and the lattice parameters of geometryoptimized structures with all functionals are very close to experimental values with a deviation <2% as shown in Table S1 in the Supporting Information. However, the variation of piezoelectric and mechanical properties with the DFT functional is small for some ZIFs, while it is large for others. Here, in the main paper, we draw only conclusions on variations in properties between different ZIFs when these variations between ZIFs are larger than the variation between functionals for the same ZIF. In the main paper, we elaborate and discuss the results obtained with B3LYP in depth as this functional has shown the best correspondence with experimental values for the elastic properties of ZIF-8. 63 Figure 3 shows the total e 14 , the clamped ion contribution (e 14 0 ), and internal strain contribution (e 14 int ) for all ZIFs. As seen in Figure Figure 3); in effect, this is a manifestation of the dipole moment of the linker changing sign going from an electrondonating to an electron-withdrawing side group. As the clamped ion contribution is a purely electronic contribution that indicates the effect of redistribution of electron density upon strain, it is informative to compare trends in the clamped ion contribution with the well-established and easily accessible parameter Hammett σ constant. It can be considered a measure of the electronic nature of the linker substituents, with a more positive Hammett constant for electron-withdrawing substituent groups, and a more negative one for the electrondonating ones. The Hammett σ constant has been published in the literature for substituted phenyl systems based on meta or para substitution as σ meta or σ para constants, respectively. 64, 65 It was also used for other chemically distinct systems such as the pK a values of imidazole derivatives and electronic, spectroscopic properties of pyrazole derivatives. 66,67 Figure 4 shows a plot of the clamped ion e 14 0 of ZIFs with Zn/Cd as the metal node and σ meta and σ para constants. Linear fit in the plot shows a good correlation with the coefficient of determination R 2 = 0.76 with σ meta and R 2 = 0.88 with σ para Hammett constants; hence, the electronic nature as captured in the Hammett constant is a good approximation to describe the effect of the linker substituent to e 14 0 in terms of changing the dipole moment of the linker. As can be seen in Figure 4, the value of e 14 0 is hardly affected by the choice of metal ion (Zn 2+ vs Cd 2+ for −CH 3 and −NO 2 substituents).

Internal Strain Contribution (e 14 int
). For all studied ZIFs, the internal strain contribution is opposite in sign and nearly equal in magnitude to the clamped ion contribution (e 14 0 ). As ZIFs are relatively soft materials, it is indeed expected that the nuclei relax (internal strain contribution) in a manner  to counter the polarization created via the clamped ion contribution. As a result, the final e 14 is much smaller for all ZIFs when compared to some inorganics such as Zn and Cd oxides and sulfides. 68,69 As shown in eq 3, the internal strain e 14 int is a product of the BECs Z s,1j * and the relaxation coefficient which depends on the change in positions of the atoms during strain. Concerning the latter, we look at changes due to external shear in all bond lengths and angles in the structure. The largest variation above a cutoff of 0.05% for the relative change in bond lengths and absolute change of 1°for angles in response to an external shear deformation of 3% is considered and shown in Table S2 in the Supporting Information. For all ZIFs, we see that the highest change in bond length occurs in the metal−N bonds ranging from 0.24% for ZIF-65 to 0.06% for ZIF-8. The variation in the bond lengths for all ZIFs being less than 0.24% indicates that bond lengths play a minor role in the shear deformation. Among the angles, the maximum variation for each ZIF constitutes that of the N−metal−N bond angle (from 2.90 to 3.48°), followed by the variation of the metal−N−C bond angle (from 0.69 to 1.44°) (see Table  S2). This indicates that the flexibility is centered around the position of the linkers relative to the Zn 2+ or Cd 2+ ions and that the organic linker is comparatively rigid. Additionally, the structural changes observed are mainly due to bond angle changes around the metal tetrahedra, while the changes in bond lengths are minimal. The other factors that influence e 14 int are the BECs, Z* which are second-rank tensors and can be represented in a matrix form in Voigt notation as In inorganics such as BaTiO 3 , anomalously high BECs are responsible for superior piezoelectric constants, where the BEC is much higher than the nominal charge of the atom (e.g., a BEC (Z Ti *) of +6.7 for Ti with a nominal charge of +4 in BaTiO 3 ). 41 Hence, we analyzed the magnitude of BECs of all atoms in the ZIF structures. For all ZIFs, the acoustic sum rule (i.e., ∑ s Z s,ij * = 0 meaning charge neutrality in all directions) is satisfied, indicating good convergence of our calculations within an accuracy of 0.07 charge units. In eq 3, BEC Z s,1j * with j = 1, 2, 3 are the relevant elements of the BEC tensor that contribute to e 14 int . Figure 5 shows a scatter plot of Z Zn/Cd,11 * , Z Zn/Cd,12 * , and Z Zn/Cd,13 * of all metal atoms Zn and Cd in the unit cell for all six ZIFs in this work. We see that Z Zn/Cd,11 * > Z Zn/Cd,12 * ,Z Zn/Cd,13 * and values of Z Zn/Cd,11 * are around +2.0 to +2.5. Thus, the overall BEC of the metal node in these ZIFs is close to their nominal charge of +2.0 and not anomalous. The influence of the linker substituent on the BEC of the metal node is negligible, with a variation of around ±0.2. The variation of the BEC of the organic atoms in discussed in the Supporting Information in Section 4.2. Overall, metal atoms have the highest BEC among all atoms in the ZIF system. In the case of the organic part of ZIFs, there is a lot of variation in the BECs depending on the displacement direction, and they mostly range from +1 to −1 indicating the covalent bond nature between them. None of the magnitudes of BECs are exceptionally high, which is one of the reasons why the internal strain component of e 14    Since the e 14 int refers to polarization along x(1) when a shear strain is applied along yz (4), this orientation to show the (011) plane was chosen in Figure 6. For the linkers aligned with their C2-axis along the x-direction, these are dipolarly aligned. Note, however, that the other linkers are oriented in a manner that diminished the overall polarization caused by the linkers.
For ZIF-8 and CdIF-1, we see a similarity in the atoms that have a higher magnitude for "A"; see Figure 6. We see that, overall, the organic linker, which as a mostly rigid unit is tilted with respect to the M 2+ tetrahedra (vide infra), contributes most to the internal strain piezoelectric coefficient e 14 int . The larger relocation of the organic atoms ( ) for these structures also, the inorganic atoms show an "A" value above the cutoff.

Elastic Compliance Constant (s 44 ) and Piezoelectric Constant (d 14 )
. The piezoelectric coefficient that determines the performance in piezoelectric devices is d 14 , which we can be derived from e 14 Figure S7 in the Supporting Information and the ratio of the standard deviation to the average varies from 0.1 to 18%, except for CdIF-1 and CdIF-8 with a deviation of 30%. Although c 44 and s 44 vary with the choice of DFT functionals for each ZIF, variation between different ZIFs is consistent in most cases. As shown in Figure S10 in Supporting Information, correspondingly there is a lot of variation of d 14 with the choice of DFT functional, due to the variability of both e 14 and s 44 with choice of functional. However, the variation between different ZIF structures (the trend) is mostly consistent. Since the B3LYP functional with and without dispersion corrections estimates reliable values close to experimentally measured mechanical properties at 295 K in ZIF-8, 63 we will discuss and compare the results of s 44 and d 14 with B3LYP for all ZIFs in the main paper. Moreover, we will only point out differences between structures when this difference is larger than the variation of the property with the functional. c 44 was previously calculated by DFT for ZIF-90 and ZIF-Cl and values are 2.502 and 3.578 GPa respectively. 70 Despite having different values for c 44 in comparison to the previous work, the values in our work computed with all the five different DFT functionals are in the same range, that is, 0.75−1.18 GPa and 0.81−0.91 GPa for ZIF-90 and ZIF-Cl, respectively (see Figure S7a).
Mechanical properties of the investigated ZIFs indicate a change in the flexibility of the framework with a metal node and linker substituent. High s 44 (or low c 44 ) means more flexibility of the framework, while low s 44 (high c 44 ) indicates stiffer mechanical properties. As seen in Figure 7, s 44 is the highest for CdIF-1 with Cd and −CH 3 as substituents and lowest for ZIF-65. Overall, we see that for the same substituent, s 44 increases significantly going from Zn 2+ to Cd 2+ and that, overall, −CH 3 is the substituent leading to the most flexibility and −NO 2 to the least. Compliance constant s 44 of ZIF-90 and ZIF-Cl range in between −CH 3 and −NO 2 substituent Zn−ZIFs. We hypothesize the cause for −NO 2 ZIFs being stiffer than −CH 3 ZIFs, although the weaker Zn−N bond in the former is due to the increased steric interactions of the larger −NO 2 groups. For example, as shown in Figure S11 in the Supporting Information, the oxygens of the −NO 2 group in adjacent linkers are closer than the hydrogens of the −CH 3 group (3.095 vs 3.297 Å where van der Waals radii of O and H are 1.5 and 1.2 Å, respectively 71 ), due to which there is more resistance during shear, resulting in higher c 44 /lower s 44 . It is clear that Cd 2+ ZIFs are more flexible with regard to the shear modulus than Zn 2+ ZIFs. This could be due to the longer Cd− N bond distance of 2.2 Å compared to the Zn−N bond distance of 2.0 Å. In the end, looking at the piezoelectric constant d 14 , it is the highest in magnitude for CdIF-1 compared to other ZIFs for all DFT functionals used in this work. Since e 14 has the same order of magnitude for all ZIFs, flexibility is the key contributing factor responsible for higher d 14 of CdIF-1. The value of d 14 for CdIF-1 is in the range of −38 to −46 pC/N (except for the M06L functional; it is −94 pC/N; see Figure  S8 in the Supporting Information). After CdIF-1, ZIF-8 seems to have d 14 values in the range of −9 to −14 pC/N. For ZIF-90, ZIF-Cl, and ZIF-65, due to huge variation in d 14 with the DFT functional, it is difficult to conclude a trend, but the magnitudes of d 14 are in the same range for these three MOFs, that is, <±8 pC/N. d 14 for CdIF-8 (marked in Figure 7) varies between 2 and 50 pC/N due to a large variation in both e 14 and s 44 with the functional, making it difficult to conclude about d 14 in CdIF-8.

Piezoelectric Properties.
Here, we compare the results of piezoelectric and mechanical properties of ZIFs in this work with some existing Zn/Cd inorganic and widely known organic piezoelectrics such as PVDF and its copolymers. Table 1 lists the compliance constants (s 44 , s pq ), clamped ion (e ik 0 ), internal strain (e ik int ), final piezoelectric constant e ik and d ik of these materials, and the ZIFs in this work.
Comparing the results of e ik 0 and e ik int for inorganic ZnO/ZnS with ZIFs, we see that inorganics have values almost an order to 2 orders of magnitude higher than ZIFs. However, one of the factors that influence e ik int is BECs and BECs of metal Z Zn/Cd * in ZIFs is similar to BECs of Zn and Cd in the inorganics presented in Table 1. This value is around +2.10 to +2.12 for Zn and + 2.13 to +2.27 for Cd in the inorganics mentioned here, and it is comparable to the values of +2.0 to +2.5 for the ZIFs. The low BECs for Zn 2+ and Cd 2+ in both ZnX and CdY ZIFs could be due to the completely occupied d orbitals, unlike Ti 4+ in BaTiO 3 with anomalous BEC which has unoccupied d orbitals. 41 To further understand the differences in e ik 0 and e ik int between ZnO/ZnS and ZIFs, we looked at the extent of clamped ion contribution that is compensated by internal strain using the ratio −e ik 0 /e ik int . For ZnO and ZnS, this is 28− 40% and 69−83%, respectively, whereas for all ZIFs, it is 80− 100%. Higher compensation in ZnS than ZnO is attributed to the bigger electronic polarizability of sulfide than the oxide anion. 40 However, for ZIFs in this work, the clamped ion contribution is much smaller than that of inorganics (∼0.06 C/ m 2 (ZIFs) versus ∼0.59 C/m 2 (ZnX)), so there is less to compensate for by the internal strain contribution. Plausible reasons for the higher compensation of clamped ion contribution are that (a) ZIFs belong to a nonpolar cubic space group and are highly symmetric and (b) their framework has higher flexibility, by which atoms relax easily in a manner to counter the polarization created via clamped ion contribution. Finally, comparing the |e 14 | of all ZIFs which is 0.01 C/m 2 while that of Zn/Cd inorganic piezoelectrics has values between 0.03 and 1.19 C/m 2 indicates a low piezoelectric constant "e". Looking at the mechanical properties, the magnitude of s 44 for all ZIFs in this study is greater than ∼560 TPa −1 and the highest s 44 among them is for CdIF-1 with a magnitude of ∼5500 TPa −1 . Between Zn/Cd inorganics and ZIFs, we see that the compliance constant s pq is very low for inorganics by almost one to 2 orders of magnitude than ZIFs. One thing commonly observed in both Zn/Cd inorganics and Zn/Cd−ZIFs is an increase in s pq or an improvement in flexibility with a change in metal from Zn to Cd.
Since data for e ik 0 and e ik int are not available in the literature for organic piezoelectrics PVDF and PVDF−trifluoroethylene (P(VDF−TrFE)), we compare here the BECs (Z*), one of the key contributors to e ik int in these materials with ZIF charges. Previous work on the ferroelectric phase β-PVDF shows that the Z* for C is −0.2/1.45 and that for F and H is −0.75 and 0.14, respectively. 73 Recent work on metal-free organic perovskite materials such as MDABCO−NH 4 −X 3 (X = Cl, Br, I) shows C, N, and H with a magnitude of Z* around 0.14 to 0.67. 11 These values are similar to BECs of the atoms in the organic linkers for ZIFs. Comparing |e 14 | of all ZIFs with PVDF and its copolymers, some e ik values are also as low as ZIFs, but few e ik values are up to an order of magnitude larger than ZIFs. Also, these organic piezoelectrics are polar materials; hence, a spontaneous polarization occurs by virtue of structure. Mechanical properties of organic materials considered here are lower than those of most ZIFs in some directions, but in a few directions, s pq is comparable with s 44 of ZIFs (except for CdIF-1 for which s 44 is 1 order of magnitude higher). The literature on clamped ion (e ik 0 ) and internal strain contributions (e ik int ) of e for Zn/Cd Se, Te, and polymer piezoelectrics considered here is not available, and these values are not included in Table 1. a Upon the overall comparison of ZIFs with Zn/Cd-based piezoelectric inorganics and specified organics, even though the piezoelectric constant e of ZIFs is lower, they are softer and flexible, which is the main contributing factor for obtaining a high piezoelectric coefficient d. This is observed in the case of CdIF-1, where |e 14 | = ∼0.01 C/m 2 , which is lower than that of most materials shown in Table 1, but it has a |d 14 | of 38 to 46 pC/N, which is higher than that of all II−VI-type inorganics and is in the same range of value as the polymer PVDF. After CdIF-1, ZIF-8 with |d 14 | = 9 to 14 pC/N is comparable to P(VDF−TrFE) and some Cd inorganics.
Comparing ZIFs with piezoelectrics other than Zn/Cd inorganics and PVDF, ZIF-8 has a d 14 comparable to d 33 of LiNbO 3 (11 pC/N), poly-L-lactic acid (PLLA) (11 pC/N), is larger than that of all the ZIFs investigated in this work, which is mainly due to larger e ik (0.35 C/m 2 ) than ZIFs because the highest s pq (650 TPa −1 ) is lower than the s 44 of ZIFs.
The study of cadmium ZIFs in this work provides a good understanding of how flexibility can influence the piezoelectric properties. However, Cd is also a heavy metal like lead, which is toxic and carcinogenic for the human body. Hence, Cd− ZIFs are a proof of concept to show the potential of MOFs as piezoelectric materials but are not ideal candidates for future applications.

Acoustic Properties.
Another factor on which the efficiency of piezoelectric devices depends is the proper matching of acoustic impedance (z) of energy harvesting materials and the media in which energy is being harvested. Hence, we discuss in brief the acoustic impedance of the ZIFs in this work and compare it with that of some piezoelectric ceramics and polymers and also the acoustic impedance of typical harvesting media. Acoustic impedance is a measure of the ease with which sound travels through a particular medium. It can be calculated as the product of acoustic velocity and density of the material. Theoretically calculated maximum longitudinal acoustic impedance values of ZIFs are shown in Table 2, and detailed maximum and minimum values of longitudinal and transverse velocities are provided in the Supporting Information in Table S3.
The acoustic impedance values of ZIFs are lower than those of both piezoelectric ceramics and organics. Among the ZIFs, CdIF-1 specifically has a lower acoustic impedance than other ZIFs. Importantly, the acoustic impedance of ZIFs is much closer to harvesting media like water and living tissues, compared to ceramics and some of the polymers. This is an advantage for harvesting energy from these media.

CONCLUSIONS
To obtain a structure−property relationship in MOFs, we systematically investigated the piezoelectric and mechanical properties of six non-centrosymmetric, nonpolar ZIFs with Zn/Cd as metal nodes and varying the substituent on the imidazolate linker (R = CH 3 , Cl, CHO, NO 2 ) using DFT calculations. Piezoelectric coefficient e is similar in magnitude for all ZIFs because of the huge compensation of clamped ion contribution e 14 0 by internal strain e 14 int , possibly from the higher flexibility of the frameworks. The magnitude of "e" for the current ZIFs ∼0.01 C/m 2 is low compared to most inorganic piezoelectrics. This can be tuned for future MOFs by (1)  Even though the e values for current ZIFs are low, the highest magnitude of piezoelectric constant d is seen for CdIF-1 with a Cd metal node and methyl (CH 3 ) imidazolate as a linker with a theoretical value of 38−46 pC/N. This value is comparable to those of the most common organic piezoelectric PVDF, and the main contributing factor for this high d value is the low shear modulus (high flexibility) of the framework. The estimated FOM of CdIF-1 calculated using an experimentally measured dielectric constant value of ZIF-8 (2.33) and d values obtained in this work is 620/ε 0 to 900/ε 0 (pC/N) 2 . This is quite high compared to the FOM ij of PZT (410 pC/N, ε = 1800) and PVDF (30 pC/N, ε = 10), 75 93/ε 0 and 90/ε 0 (pC/N) 2 , respectively. Hence, when tunable parameters such as structural building units are optimized to obtain high flexibility, high BECs, and favorable non-centrosymmetric symmetry, MOFs should be exceptional candidates to be used as piezoelectric energy harvesters.